# Call the data cleaning script
source("./code/data_cleaning.R")Attached: 'Groundhog' (Version: 3.2.0)
Tips and troubleshooting: https://groundhogR.com
The following tracks the data analysis for our 2024 Undergraduate Student Research Award (USRA) project, “What accounts for the relationship between disgust and religiosity?” This is the master analysis document, so all analyses are documented here—with any other code called within a chunk below. The goals of this analysis are to (1) clean the data, (2) check for internal reliability of our measures, (3) generate descriptive statistics, (4) visualize the distribution of our variables, (5) test our hypotheses, and (6) conduct exploratory analyses as appropriate.
Individual differences in disgust sensitivity have repeatedly been shown to relate to religiosity, as predicted by the behavioral immune system model of disgust, leading to some scholars suggesting that religiosity is an evolved disease avoidance strategy. However, there are multiple mechanisms that aim to account for how religion contributes to disease avoidance, and these different mechanisms provide predictions for the psychological mediators of the relationship between disgust and religiosity. This study aims to test three such mediators, to better differentiate and understand why disgust sensitivity relates to religiosity. These mediators are represented by the following hypotheses.
Strong version of statistical hypothesis:
Weak version of statistical hypothesis:
Strong version of statistical hypothesis:
Weak version of statistical hypothesis:
Strong version of statistical hypothesis:
Weak version of statistical hypothesis:
We will test these hypotheses through mediation analysis, the specifics of which will be reviewed in the Testing Hypotheses section.
All data files are saved in .csv format in the data directory, and the steps that were taken to process these files are documented in data/data_preparation_notes.txt. To clean the data for analysis, we will call a pre-written script (code/data_cleaning.R) that returns a data frame called data. The script also saves this data frame as data/data_clean.csv.
# Call the data cleaning script
source("./code/data_cleaning.R")Attached: 'Groundhog' (Version: 3.2.0)
Tips and troubleshooting: https://groundhogR.com
Note: The Three Domains of Disgust Scale (TDDS) is generally administered with seven response options (Olatunji et al., 2012), only anchored with 0 (Not disgusting at all) and 6 (Extremely disgusting). Due to a survey coding error, our sample completed the TDDS with only six response options and the same anchors. Because of this, each TDDS item has a range from 0 to 5 (rather than 6), and raw scores for the full scale and subscales will be underestimated accordingly.
The data data frame contains the following variables (and column names).
Participant ID (ID):
Start Time (start_time) and Finish Time (finish_time):
Age (age):
Sex (sex):
Gender (gender):
The gender of the participant. Responses were standardized into six categories formulated based on the genders reported. To see which original responses were put into these bins, see the appropriate section of the code/data_cleaning.R script. The categories are as follows:
Ethnicity (ethnicity):
The ethnicity of the participant. Participants were provided with a list of potential ethnicities to choose from as well as the options of writing in their own or “Rather Not Say”. Those who reported multiple ethnicities were categorized as “Mixed”. There are 11 unique ethnicities for this sample:
Nationality (nationality):
Religious Affiliation (religious_affiliation):
Christian Affiliation (christian_affiliation):
Religious Importance (religious_importance):
The Centrality of Religiosity Scale (CRS):
The Three Domains of Disgust Scale (full, TDDS_f; pathogen, TDDS_p; sexual, TDDS_s, moral, TDDS_m):
The TDDS was developed by Tybur et al. (2009) in response to a lack of theoretical and empirical grounding for prior measures of disgust sensitivity. Based on evolutionary theory, they predicted that items which people find disgusting should form three factors, one for pathogen, sexual, and moral disgust. Across multiple studies, they demonstrate that a three-factor solution is most parsimonious and fits the data well, that the full scale has good concurrent validity, and that the subscales have good convergent and divergent validity. Subsequent investigations have confirmed the three factor structure of the scale, although they find poorer validity evidence for the moral disgust subscale (Olatunji et al., 2012). All three subscales will be used in exploratory analyses and bivariate correlations, but the pathogen disgust scale will be used to test our central three hypotheses and six specific statistical hypotheses.
Scores for the full scale and each subscale were calculated by summing item values.
The Revised Sociosexual Orientation Inventory (full, SOI_f; behavior, SOI_b; attitudes, SOI_a; desire, SOI_d):
The original Sociosexual Orientation Inventory (SOI) was a 7-item global measure of sociosexual orientation (Simpson & Gangestad, 1991), but this measure included behavioral, attitudinal, and (one) desire items together, with different response formats. This is problematic because each of these three components are differentiable both conceptually and empirically, and the different response formats often lead to an improper weighting of each of the three components contributing to the total score, resulting in poor reliability. Penke & Asendorpf (2008) developed the SOI-R to address these, and other, problems. The SOI-R is a 9-item measure with three items each tailored towards behaviors, attitudes, and desires, which they differentiated based on exploratory and confirmatory factor analysis as well as convergent and divergent validation. Each of the subscales show adequate reliability, along with the global score. As we are interested in operationalizing a monogamous mating strategy, it is participants intentions that we were most interested in measuring. In this case, one’s desire and behavior may, or may not, correspond to their intentions, which would best be operationalized as sociosexual attitudes. Therefore, we will use the attitudes subscale of the SOI-R to operationalize a monogamous mating strategy. This will be used to test our mediational hypothesis involving monogamous sexual strategies.
Scores for the full scale and each subscale were calculated by taking the arithmetic mean of item values.
The In-Group Preference Subscale of the Short Form of the Generalized Ethnocentrism Scale (SFGENE-7; full, GENE_f; preference, GENE_p; superiority, GENE_s):
The SFGENE-7 is simply a short form of the Generalized Ethnocentrism Scale (Neuliep, 2002; Neuliep & McCroskey, 1997). The GENE is composed of 22 items, 15 of which are scored, and out of all the measures that we have, the GENE would be by far the longest. By reducing the number of items to 7 with the SFGENE-7, we can allow for the recruitment of more participants. The SFGENE-7 has shown adequate reliability despite a loss of items, as well as evidence for convergent validity with constructs like tolerance and multicultural ideology. In addition, it has been broken down via factor analysis into two distinct subscales, the first three items comprising the in-group preference scale and the last four items comprosing the in-group superiority scale. Each of these subscales also have strong internal reliabilities, and the in-group preference scale shows stronger negative correlations with tolerance and multicultural identity than the in-group superiority scale. We will use use only the first three items to measure the in-group preference aspect of ethnocentrism in order to operationalize out-group avoidance, particularly because the items comprising the in-group superiority subscale are quite similar at face value to items in the conventionalism scale that we are using to measure the tendency towards adherence to tradition. The in-group preference subscale will be used to test the mediational hypothesis involving out-group avoidance.
Items were summed to provide the total, preference, and superiority scores.
Conventionalism Subscale of Aggression-Submission-Conventionalism Scale of Authoritarianism (ASC-C; CONV) and additional Traditionalism Items (TRAD; together, CONV_f):
TRAD items in the data frame. We will calculate Chronbach’s alpha and only accept these items in the measure if we keep an alpha above .8. If the alpha is below .8, we will remove items until it performs at this level.To quickly assess the internal reliability of our scales, we will calculate Chronbach’s alpha for each of our scales. We will start with the CONV_f scale, because we may remove items we added to it.
# Defining the packages to use to do the alpha analysis
pkg <- c("tidyverse", "psych")
# tidyverse if not attached from data cleaning script (for pipe and select function)
# psych for alpha (and later descriptives)
# Loading the groundhog package to install packages from a certain date
library(groundhog)
# Reading in the packages with the groundhog package for reproducibility
# Message suppressed in order to allow for rendering of quarto document
suppressMessages(groundhog.library(pkg = pkg, date = "2024-06-01"))
# Remove the pkg vector
rm(pkg)# Chronbach's Alpha for CONV_f (which includes all CONV and TRAD items)
alpha(x = data %>% select(starts_with("CONV."), starts_with("TRAD.")))
Reliability analysis
Call: alpha(x = data %>% select(starts_with("CONV."), starts_with("TRAD.")))
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.83 0.83 0.87 0.29 4.9 0.015 3 0.86 0.28
95% confidence boundaries
lower alpha upper
Feldt 0.8 0.83 0.86
Duhachek 0.8 0.83 0.86
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
CONV.1 0.82 0.82 0.86 0.30 4.6 0.015 0.034 0.29
CONV.2 0.80 0.80 0.84 0.27 4.0 0.017 0.031 0.25
CONV.3 0.82 0.82 0.85 0.29 4.6 0.015 0.029 0.30
CONV.4 0.81 0.81 0.85 0.28 4.4 0.016 0.034 0.28
CONV.5 0.82 0.81 0.85 0.28 4.4 0.016 0.034 0.28
CONV.6 0.82 0.82 0.85 0.29 4.5 0.016 0.028 0.29
TRAD.1 0.84 0.84 0.87 0.33 5.4 0.014 0.023 0.31
TRAD.2 0.82 0.82 0.85 0.29 4.5 0.016 0.030 0.28
TRAD.3 0.82 0.82 0.86 0.29 4.5 0.015 0.034 0.28
TRAD.4 0.82 0.81 0.85 0.28 4.4 0.016 0.031 0.29
TRAD.5 0.81 0.80 0.84 0.27 4.1 0.017 0.028 0.26
TRAD.6 0.81 0.81 0.85 0.28 4.3 0.016 0.030 0.28
Item statistics
n raw.r std.r r.cor r.drop mean sd
CONV.1 289 0.53 0.53 0.47 0.42 2.3 1.4
CONV.2 289 0.76 0.75 0.74 0.68 3.4 1.5
CONV.3 289 0.54 0.55 0.52 0.43 3.4 1.5
CONV.4 289 0.64 0.63 0.58 0.53 2.7 1.5
CONV.5 289 0.63 0.63 0.59 0.54 2.2 1.5
CONV.6 289 0.57 0.58 0.55 0.47 3.9 1.3
TRAD.1 289 0.27 0.27 0.17 0.14 1.5 1.4
TRAD.2 289 0.57 0.57 0.53 0.47 3.9 1.4
TRAD.3 289 0.57 0.57 0.51 0.46 2.0 1.5
TRAD.4 289 0.63 0.63 0.60 0.53 3.0 1.5
TRAD.5 289 0.73 0.73 0.73 0.66 3.7 1.4
TRAD.6 289 0.63 0.64 0.61 0.54 3.7 1.4
Non missing response frequency for each item
0 1 2 3 4 5 6 miss
CONV.1 0.13 0.17 0.27 0.25 0.12 0.03 0.02 0
CONV.2 0.04 0.07 0.14 0.28 0.24 0.13 0.10 0
CONV.3 0.04 0.06 0.12 0.29 0.25 0.15 0.09 0
CONV.4 0.07 0.13 0.27 0.26 0.13 0.08 0.06 0
CONV.5 0.13 0.18 0.31 0.21 0.10 0.04 0.03 0
CONV.6 0.01 0.03 0.09 0.25 0.30 0.19 0.12 0
TRAD.1 0.31 0.25 0.22 0.15 0.03 0.02 0.01 0
TRAD.2 0.02 0.04 0.10 0.20 0.31 0.17 0.15 0
TRAD.3 0.19 0.23 0.25 0.16 0.11 0.04 0.01 0
TRAD.4 0.06 0.09 0.23 0.25 0.23 0.09 0.06 0
TRAD.5 0.03 0.05 0.09 0.26 0.32 0.13 0.12 0
TRAD.6 0.02 0.07 0.08 0.25 0.30 0.20 0.08 0
The raw alpha value (.831; standardized .830) is acceptable here, so we will use our full CONV_f variable to operationalize traditionalism.
To compare these values to those for the original conventionalism measure, I will calculate its’ alpha as well.
# Chronbach's Alpha for CONV (only CONV items)
alpha(x = data %>% select(starts_with("CONV.")))
Reliability analysis
Call: alpha(x = data %>% select(starts_with("CONV.")))
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.74 0.74 0.76 0.33 2.9 0.024 3 0.97 0.38
95% confidence boundaries
lower alpha upper
Feldt 0.7 0.74 0.79
Duhachek 0.7 0.74 0.79
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
CONV.1 0.73 0.73 0.75 0.36 2.8 0.025 0.028 0.40
CONV.2 0.67 0.67 0.70 0.29 2.1 0.031 0.038 0.21
CONV.3 0.72 0.72 0.70 0.34 2.6 0.025 0.018 0.39
CONV.4 0.70 0.70 0.71 0.32 2.4 0.028 0.031 0.36
CONV.5 0.69 0.69 0.70 0.31 2.3 0.029 0.030 0.35
CONV.6 0.72 0.71 0.70 0.33 2.5 0.026 0.020 0.39
Item statistics
n raw.r std.r r.cor r.drop mean sd
CONV.1 289 0.59 0.59 0.45 0.39 2.3 1.4
CONV.2 289 0.75 0.75 0.67 0.59 3.4 1.5
CONV.3 289 0.62 0.63 0.56 0.43 3.4 1.5
CONV.4 289 0.68 0.67 0.58 0.50 2.7 1.5
CONV.5 289 0.70 0.70 0.62 0.53 2.2 1.5
CONV.6 289 0.63 0.64 0.58 0.45 3.9 1.3
Non missing response frequency for each item
0 1 2 3 4 5 6 miss
CONV.1 0.13 0.17 0.27 0.25 0.12 0.03 0.02 0
CONV.2 0.04 0.07 0.14 0.28 0.24 0.13 0.10 0
CONV.3 0.04 0.06 0.12 0.29 0.25 0.15 0.09 0
CONV.4 0.07 0.13 0.27 0.26 0.13 0.08 0.06 0
CONV.5 0.13 0.18 0.31 0.21 0.10 0.04 0.03 0
CONV.6 0.01 0.03 0.09 0.25 0.30 0.19 0.12 0
CONV_f variable has stronger internal reliability than the CONV items alone.# Chronbach's Alpha for Centrality of Religiosity Items
alpha(x = data %>% select(starts_with("CRS.")))
Reliability analysis
Call: alpha(x = data %>% select(starts_with("CRS.")))
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.92 0.92 0.91 0.69 11 0.0069 2.6 1.3 0.72
95% confidence boundaries
lower alpha upper
Feldt 0.90 0.92 0.93
Duhachek 0.91 0.92 0.93
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
CRS.1 0.93 0.93 0.92 0.78 13.8 0.0064 0.002 0.80
CRS.2 0.89 0.89 0.88 0.68 8.4 0.0094 0.012 0.66
CRS.3 0.90 0.90 0.88 0.69 8.8 0.0090 0.016 0.70
CRS.4 0.89 0.89 0.86 0.66 7.7 0.0104 0.011 0.65
CRS.5 0.89 0.89 0.87 0.67 8.1 0.0096 0.015 0.65
Item statistics
n raw.r std.r r.cor r.drop mean sd
CRS.1 289 0.73 0.76 0.65 0.63 2.6 1.1
CRS.2 289 0.90 0.89 0.87 0.83 3.0 1.5
CRS.3 289 0.88 0.88 0.84 0.81 2.3 1.4
CRS.4 289 0.93 0.92 0.91 0.87 2.5 1.7
CRS.5 289 0.90 0.90 0.88 0.85 2.5 1.4
Non missing response frequency for each item
1 2 3 4 5 miss
CRS.1 0.14 0.37 0.28 0.13 0.08 0
CRS.2 0.18 0.29 0.15 0.11 0.27 0
CRS.3 0.42 0.22 0.13 0.09 0.14 0
CRS.4 0.46 0.15 0.08 0.03 0.28 0
CRS.5 0.33 0.26 0.16 0.11 0.14 0
# Chronbach's Alpha for Pathogen Disgust Items
alpha(x = data %>% select(., c("TDDS.12", "TDDS.15", "TDDS.9", "TDDS.3", "TDDS.21", "TDDS.18", "TDDS.6")))
Reliability analysis
Call: alpha(x = data %>% select(., c("TDDS.12", "TDDS.15", "TDDS.9",
"TDDS.3", "TDDS.21", "TDDS.18", "TDDS.6")))
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.77 0.77 0.76 0.33 3.4 0.021 3 0.91 0.34
95% confidence boundaries
lower alpha upper
Feldt 0.73 0.77 0.81
Duhachek 0.73 0.77 0.81
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
TDDS.12 0.74 0.74 0.72 0.33 2.9 0.023 0.0094 0.34
TDDS.15 0.72 0.72 0.70 0.30 2.6 0.026 0.0063 0.28
TDDS.9 0.72 0.73 0.70 0.31 2.7 0.025 0.0064 0.31
TDDS.3 0.75 0.75 0.74 0.34 3.1 0.023 0.0091 0.34
TDDS.21 0.75 0.75 0.74 0.33 3.0 0.023 0.0114 0.35
TDDS.18 0.75 0.75 0.74 0.34 3.1 0.022 0.0095 0.35
TDDS.6 0.75 0.76 0.74 0.34 3.1 0.023 0.0076 0.34
Item statistics
n raw.r std.r r.cor r.drop mean sd
TDDS.12 289 0.65 0.65 0.57 0.49 3.4 1.4
TDDS.15 289 0.73 0.74 0.70 0.61 3.3 1.3
TDDS.9 289 0.71 0.71 0.67 0.57 2.5 1.4
TDDS.3 289 0.58 0.61 0.51 0.44 4.0 1.1
TDDS.21 289 0.65 0.63 0.52 0.47 3.0 1.5
TDDS.18 289 0.63 0.61 0.51 0.44 3.1 1.5
TDDS.6 289 0.61 0.60 0.50 0.44 1.9 1.4
Non missing response frequency for each item
0 1 2 3 4 5 miss
TDDS.12 0.02 0.09 0.15 0.21 0.25 0.28 0
TDDS.15 0.01 0.11 0.15 0.26 0.24 0.24 0
TDDS.9 0.07 0.20 0.28 0.21 0.15 0.10 0
TDDS.3 0.00 0.03 0.09 0.17 0.23 0.48 0
TDDS.21 0.06 0.15 0.18 0.22 0.15 0.24 0
TDDS.18 0.04 0.14 0.18 0.19 0.17 0.27 0
TDDS.6 0.17 0.24 0.26 0.18 0.08 0.06 0
TDDS_p) shows adequate internal reliability here (alpha = .769; standardized = .772).# Chronbach's Alpha for the Attitudes Subscale of the SOI-R
alpha(x = data %>% select(., c("SOI.4", "SOI.5", "SOI.6")))
Reliability analysis
Call: alpha(x = data %>% select(., c("SOI.4", "SOI.5", "SOI.6")))
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.81 0.81 0.74 0.58 4.2 0.02 5.2 2.4 0.6
95% confidence boundaries
lower alpha upper
Feldt 0.76 0.81 0.84
Duhachek 0.77 0.81 0.84
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
SOI.4 0.76 0.76 0.62 0.62 3.2 0.028 NA 0.62
SOI.5 0.69 0.69 0.52 0.52 2.2 0.037 NA 0.52
SOI.6 0.75 0.75 0.60 0.60 3.0 0.029 NA 0.60
Item statistics
n raw.r std.r r.cor r.drop mean sd
SOI.4 289 0.83 0.83 0.70 0.63 5.9 2.7
SOI.5 289 0.88 0.87 0.78 0.70 4.3 2.9
SOI.6 289 0.84 0.84 0.71 0.64 5.3 2.8
Non missing response frequency for each item
1 2 3 4 5 6 7 8 9 miss
SOI.4 0.11 0.04 0.08 0.06 0.11 0.09 0.15 0.11 0.25 0
SOI.5 0.29 0.11 0.10 0.04 0.10 0.09 0.09 0.06 0.13 0
SOI.6 0.17 0.04 0.11 0.05 0.14 0.09 0.08 0.12 0.19 0
SOI_a) seems to have good internal consistency here (alpha = .806; standardized = .806).# Chronbach's Alpha for the Preferences Subscale of the GENE
alpha(x = data %>% select(., c("GENE.1", "GENE.2", "GENE.3")))
Reliability analysis
Call: alpha(x = data %>% select(., c("GENE.1", "GENE.2", "GENE.3")))
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.52 0.53 0.55 0.27 1.1 0.05 2.3 0.68 0.095
95% confidence boundaries
lower alpha upper
Feldt 0.42 0.52 0.61
Duhachek 0.43 0.52 0.62
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
GENE.1 0.801 0.804 0.672 0.672 4.09 0.023 NA 0.672
GENE.2 0.093 0.093 0.049 0.049 0.10 0.107 NA 0.049
GENE.3 0.174 0.174 0.095 0.095 0.21 0.097 NA 0.095
Item statistics
n raw.r std.r r.cor r.drop mean sd
GENE.1 289 0.54 0.53 0.096 0.078 3.2 0.95
GENE.2 289 0.81 0.82 0.766 0.536 1.7 0.89
GENE.3 289 0.80 0.80 0.734 0.472 1.9 0.99
Non missing response frequency for each item
1 2 3 4 5 miss
GENE.1 0.05 0.12 0.47 0.27 0.09 0
GENE.2 0.51 0.32 0.12 0.03 0.01 0
GENE.3 0.47 0.29 0.18 0.06 0.01 0
GENE_p_2.# Creating the GENE_p_2 variable as the sum of items 2 and 3
data <- data %>%
mutate(
GENE_p_2 = rowSums(select(., c(GENE.2, GENE.3)))
)First, we will see how many participants we have and how long it took participants to complete the survey.
# Calculate the number of participants in the sample
nrow(data)[1] 289
# Calculate descriptive statistics for the difference between the finish and start time (in minutes)
describe(as.numeric(difftime(data$finish_time, data$start_time, units = "mins"))) vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 289 14.39 58.1 8 8.94 2.97 3 821 818 12.25 154.08 3.42
In total, we have N = 289 participants with complete data.
Although the mean time taken is around 14.39, the SD is huge, and the maximum is a massive 821 minutes (over 13 hours). The trimmed mean of 8.94 is therefore a better estimate.
Next, we will describe the sample in terms of the categorical variables. The following generates frequency tables for all categorical variables.
# Create a frequency table for each factor variable
lapply(data[sapply(data, is.factor)], table)$sex
Male Female
146 143
$gender
Agender Female Gender Fluid Gender Queer Male Non-binary
2 135 1 1 144 4
$ethnicity
African Black or African American
52 4
Caribbean East Asian
0 14
Latino or Hispanic Middle Eastern
62 4
Mixed Native American or Alaskan Native
11 0
South Asian White or Causasian
8 122
White or Sapharic Jew Black British
6 1
White Mexican Romani or Traveller
0 0
South East Asian Rather Not Say
2 3
$nationality
Algeria Australia Austria Belgium Brazil
1 7 1 1 3
Canada Chile China Columbia Czech Republic
30 13 1 2 5
Egypt Eritrea Finland France Germany
1 1 2 3 4
Ghana Greece Hungary India Ireland
1 5 3 1 1
Israel Italy Japan Kenya Mexico
2 5 1 6 39
New Zealand Poland Portugal Romania Saudi Arabia
6 14 25 1 1
Slovakia South Africa Spain Sweden Tunisia
1 46 7 5 1
Turkey United Kingdom United States Venezuela Vietnam
2 21 17 1 1
Zimbabwe
1
$religious_affiliation
Agnostic Anti-religious Atheist Buddhist Christian
6 1 2 2 123
Deist Hindu Muslim None Pagan
1 1 11 139 1
Spiritual Spiritualist
1 1
$christian_affiliation
Anglican Apostolic Baptist
1 1 11
Catholic Congregationalist Evangelical
68 1 4
Jehovah’s Witness Lutheran Methodist
1 6 9
Non-religious Christian None Orthodox
1 4 4
Pentecostal Presbyterian Protestant
1 2 8
Seventh Day Adventist
1
The output is likely difficult to read, so I will break it down here.
Sex:
Gender:
Male: 144 (49.8%)
Female: 135 (46.7%)
Non-Binary: 4 (1.4%)
Agender: 2 (.7%)
Gender Fluid: 1 (.3%)
Gender Queer: 1 (.3%)
Ethnicity:
Nationality:
Religious Affiliation:
Christian Affiliation:
Now we will describe our numeric variables starting with our only demographic numeric variable, age.
# Generate descriptive statistics for age
describe(data$age) vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 289 30.85 10.85 28 29.18 7.41 18 86 68 1.74 3.84 0.64
Now we will describe the important scales for our analysis.
# Describe important scales in our analysis
describe(data %>% select("religious_importance", "CRS", "TDDS_f", "TDDS_p", "SOI_f", "SOI_a", "GENE_p", "GENE_p_2", "CONV_f")) vars n mean sd median trimmed mad min max range
religious_importance 1 289 30.70 35.63 12.00 26.08 17.79 0 100 100
CRS 2 289 2.59 1.26 2.20 2.50 1.19 1 5 4
TDDS_f 3 289 60.12 16.05 59.00 59.82 16.31 12 100 88
TDDS_p 4 289 21.26 6.35 22.00 21.30 5.93 4 35 31
SOI_f 5 289 3.58 1.62 3.22 3.48 1.81 1 8 7
SOI_a 6 289 5.16 2.40 5.00 5.19 2.47 1 9 8
GENE_p 7 289 6.81 2.03 7.00 6.72 1.48 3 15 12
GENE_p_2 8 289 3.57 1.71 3.00 3.34 1.48 2 10 8
CONV_f 9 289 35.55 10.30 36.00 35.78 10.38 0 65 65
skew kurtosis se
religious_importance 0.87 -0.83 2.10
CRS 0.57 -1.06 0.07
TDDS_f 0.17 -0.29 0.94
TDDS_p -0.07 -0.49 0.37
SOI_f 0.53 -0.59 0.10
SOI_a -0.03 -1.07 0.14
GENE_p 0.50 0.36 0.12
GENE_p_2 1.02 0.64 0.10
CONV_f -0.24 0.44 0.61
religious_importance and the GENE_p_2, there does not seem to be a problem here. Variables will be inspected more closely when testing assumptions for statistical tests.Before moving on, we will calculate correlations between each of these variables.
# Calculate correlations for numerical variables
corr.test(data %>% select("age", "religious_importance", "CRS", "TDDS_f", "TDDS_p", "TDDS_s", "SOI_f", "SOI_a", "GENE_p", "GENE_p_2", "CONV_f"))Call:corr.test(x = data %>% select("age", "religious_importance",
"CRS", "TDDS_f", "TDDS_p", "TDDS_s", "SOI_f", "SOI_a", "GENE_p",
"GENE_p_2", "CONV_f"))
Correlation matrix
age religious_importance CRS TDDS_f TDDS_p TDDS_s
age 1.00 -0.07 -0.09 0.02 -0.08 -0.09
religious_importance -0.07 1.00 0.92 0.46 0.22 0.52
CRS -0.09 0.92 1.00 0.47 0.24 0.53
TDDS_f 0.02 0.46 0.47 1.00 0.76 0.82
TDDS_p -0.08 0.22 0.24 0.76 1.00 0.47
TDDS_s -0.09 0.52 0.53 0.82 0.47 1.00
SOI_f 0.12 -0.28 -0.27 -0.38 -0.16 -0.54
SOI_a 0.14 -0.39 -0.41 -0.45 -0.22 -0.59
GENE_p -0.06 0.07 0.07 0.02 0.01 0.07
GENE_p_2 0.03 0.02 0.00 -0.02 -0.03 0.04
CONV_f 0.08 0.33 0.31 0.12 0.01 0.13
SOI_f SOI_a GENE_p GENE_p_2 CONV_f
age 0.12 0.14 -0.06 0.03 0.08
religious_importance -0.28 -0.39 0.07 0.02 0.33
CRS -0.27 -0.41 0.07 0.00 0.31
TDDS_f -0.38 -0.45 0.02 -0.02 0.12
TDDS_p -0.16 -0.22 0.01 -0.03 0.01
TDDS_s -0.54 -0.59 0.07 0.04 0.13
SOI_f 1.00 0.88 -0.11 -0.10 -0.10
SOI_a 0.88 1.00 -0.15 -0.11 -0.15
GENE_p -0.11 -0.15 1.00 0.88 0.08
GENE_p_2 -0.10 -0.11 0.88 1.00 0.03
CONV_f -0.10 -0.15 0.08 0.03 1.00
Sample Size
[1] 289
Probability values (Entries above the diagonal are adjusted for multiple tests.)
age religious_importance CRS TDDS_f TDDS_p TDDS_s SOI_f
age 0.00 1.00 1.00 1.00 1.00 1.00 1.00
religious_importance 0.21 0.00 0.00 0.00 0.01 0.00 0.00
CRS 0.15 0.00 0.00 0.00 0.00 0.00 0.00
TDDS_f 0.79 0.00 0.00 0.00 0.00 0.00 0.00
TDDS_p 0.15 0.00 0.00 0.00 0.00 0.00 0.21
TDDS_s 0.14 0.00 0.00 0.00 0.00 0.00 0.00
SOI_f 0.05 0.00 0.00 0.00 0.01 0.00 0.00
SOI_a 0.01 0.00 0.00 0.00 0.00 0.00 0.00
GENE_p 0.35 0.20 0.25 0.75 0.81 0.21 0.06
GENE_p_2 0.66 0.79 1.00 0.72 0.67 0.55 0.09
CONV_f 0.17 0.00 0.00 0.04 0.93 0.03 0.10
SOI_a GENE_p GENE_p_2 CONV_f
age 0.42 1.00 1.00 1.00
religious_importance 0.00 1.00 1.00 0.00
CRS 0.00 1.00 1.00 0.00
TDDS_f 0.00 1.00 1.00 1.00
TDDS_p 0.01 1.00 1.00 1.00
TDDS_s 0.00 1.00 1.00 0.74
SOI_f 0.00 1.00 1.00 1.00
SOI_a 0.00 0.29 1.00 0.32
GENE_p 0.01 0.00 0.00 1.00
GENE_p_2 0.06 0.00 0.00 1.00
CONV_f 0.01 0.16 0.62 0.00
To see confidence intervals of the correlations, print with the short=FALSE option
For each correlation the sample size is N = 289.
Noteworthy for this analysis, pathogen disgust is positively and significantly related to both religious_importance (r = .22, p < .01) and the CRS (r = .24, p < .01), but sexual disgust is more strongly related to both religious_importance (r = .52, p < .01) and the CRS (r = .53, p < .01). This replicates previous work and provides initial support for the sexual strategies hypothesis.
In addition, there is a very strong correlation between religious_importance and the CRS (r = .92, p < .01).
However, neither the GENE_p nor the GENE_p_2 are significantly related to either religious_importance (r = .07, p = .20; r = .02, p = .02; respectively) or the CRS (r = .07, p = .25; r = .00, p = .00; respectively), suggesting that it could not mediate the relationship between disgust sensitivity and religiosity.
Unlike the GENE variables, CONV_f is significantly related to both religious_importance (r = .33, p < .01) and the CRS (r = .31, p < .01).
Next, the only potential mediator that is significantly related to pathogen disgust is the SOI_a:
SOI_a: r = -.22, p = .01
GENE_p: r = .01, p = 1
GENE_p_2: r = -.03, p = 1CONV_f: r = .01, p = 1
This suggests that the only plausible mediator (here, that is) of the relationship between disgust and religiosity could be monogamous mating strategies.
Now we will visualization the distribution of our variables by creating bar plots and histograms (for categorical and numerical variables, respectively).
# Loop through each categorical variable in the data frame and create a bar chart
for (var in names(data[, sapply(data, is.factor)])) {
p <- ggplot(data[, sapply(data, is.factor)], aes(x = !!sym(var))) +
geom_bar() +
theme_minimal() +
labs(title = paste("Bar Chart of", var), x = var, y = "Count") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 8)) # Adjust x axis labels to make them readable (not all over eachother)
# Print each bar chart
print(p)
# Remove the plot object so it doesn't clutter the environment
rm(p)
}# Remove the var object
rm(var)# Loop through numeric variable identified to produce a histogram
for (var in names(data %>% select("age", "religious_importance", "CRS", "TDDS_f", "TDDS_p", "SOI_f", "SOI_a", "GENE_p", "GENE_p_2", "CONV_f"))) {
p <- ggplot(data %>% select("age", "religious_importance", "CRS", "TDDS_f", "TDDS_p", "SOI_f", "SOI_a", "GENE_p", "GENE_p_2", "CONV_f"), aes(x = .data[[var]])) +
geom_histogram(bins = 10, fill = "steelblue", color = "black") +
theme_minimal() +
labs(title = paste("Histogram of", var), x = var, y = "Frequency")
# Print each histogram
print(p)
# Remove the p object so it doesn't clutter up the environment
rm(p)
}# Remove the var object from the environment
rm(var)As above, the following are our hypotheses.
Strong version of statistical hypothesis:
Weak version of statistical hypothesis:
Strong version of statistical hypothesis:
Weak version of statistical hypothesis:
Strong version of statistical hypothesis:
Weak version of statistical hypothesis:
The three hypotheses are each broken down into two statistical hypotheses, which represent their strong and weak version. In the strong version there is full mediation, and in the weak version there is only partial mediation. We have defined and preregistered inference criteria for both of these versions of each hypothesis.
Strong Version Inference Criterion:
For each model (each hypothesis) we will use the “K” method outlined by Beribisky et al. (2020) to infer full mediation. That is, if the proportion of variance explained by the indirect effect compared to the total effect is 80% or higher, we will infer full mediation. This will provide support for the strong version of the respective hypothesis.
Weak Version Inference Criterion:
Where there is not full mediation, for each model (each hypothesis) if the 95% confidence interval for the indirect effect of parasite disgust sensitivity on religiosity through the mediator does not contain zero, then we will infer partial mediation. This will provide support for the weak version of the respective hypothesis.
Before testing our hypotheses, we will standardize each variable. We will also reverse the scores of SOI_a before standardizing it, because higher scores currently represent more unrestricted sociosexual attitudes, whereas we would like higher scores to represent more restricted sociosexual attitudes (as a proxy for a monogamous mating strategy).
# Reversing the scoring for the SOI_a (saving it as SOI_a_r) and standardizing all variables
data <- data %>%
mutate(
SOI_a_r = (max(data$SOI_a) + 1) - SOI_a,
SOI_a_r_z = scale(SOI_a_r),
CRS_z = scale(CRS),
religious_importance_z = scale(religious_importance),
GENE_p_z = scale(GENE_p),
GENE_p_2_z = scale(GENE_p_2),
CONV_f_z = scale(CONV_f),
TDDS_p_z = scale(TDDS_p)
)
# Ensure that the variables are plain numeric vectors
data$TDDS_p_z <- as.numeric(data$TDDS_p_z)
data$CONV_f_z <- as.numeric(data$CONV_f_z)
data$CRS_z <- as.numeric(data$CRS_z)
data$SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
data$religious_importance_z <- as.numeric(data$religious_importance_z)
data$GENE_p_z <- as.numeric(data$GENE_p_z)
data$GENE_p_2_z <- as.numeric(data$GENE_p_2_z)
# Because the boot function was throwing a fit about this earlier
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)We will also load the packages involved in the mediation analysis.
# Defining the packages to use to do the mediation analysis
pkg <- c("car", "mediation")
# car for testing for multicollinearity and bootstrapping
# mediation for testing the indirect effect with bootstrapped CIs and accessing
# Reading in the packages with the groundhog package for reproducibility
# Message suppressed in order to allow for rendering of quarto document
suppressMessages(groundhog.library(pkg = pkg, date = "2024-06-01"))
# Dropping the pkg vector
rm(pkg)Before moving on, we will first test the baseline model where we regress the dependent variable (CRS) on the independent variable (TDDS_p), because each mediation analysis requires that there is a significant positive relationship between these variables. We will also test each of the assumptions.
# Fitting the model
baseline_model <- lm(CRS_z ~ TDDS_p_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(baseline_model, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(baseline_model))# QQ-plot for normality of residuals
qqnorm(residuals(baseline_model))
qqline(residuals(baseline_model), col = "red")The histogram of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just around -1 SD and around 2 SD. The two clusters of similar errors in prediction may be due to sex differences in disgust and religiosity. In exploratory analyses, we will add sex as a control variable.
Because our inferential criteria do not require parametric tests, we should not be in too much trouble with a non-normality of residuals. We were relying on a parametric t-test to assess whether pathogen disgust is related to religiosity in the first place, however. Because we have failed the assumption of normality of residuals, we will use a bootstrapped test to make our inference here as well.
Multicollinearity: Not applicable becuase only one predictor
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(baseline_model, col = "black", lwd = 1)# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_baseline_model <- Boot(baseline_model, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_baseline_model) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 1.6683e-17 -0.00066792 0.058534 -0.001485
TDDS_p_z 2.4299e-01 0.00007957 0.059718 0.245364
confint(boot_baseline_model) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1064481 0.1239702
TDDS_p_z 0.1206344 0.3523252
# Removing the CRS_z and TDDS_p_z vectors in the environment that the Boot package required outside of the data.frame
rm(CRS_z, TDDS_p_z)We can see from the Scatter Plot that there is a positive relationship between pathogen disgust and religiosity.
Based on the estimates, the standardized coefficient for pathogen disgust is around .24. This aligns with estimates from a meta-analysis of this relationship (Yu et al., 2022).
This shows that the IV effects the DV, as in the second regression equation for testing mediation in Baron & Kenny (1986). To establish mediation, we will still need to establish that the IV affects the mediator and the mediator affects the DV. Before we test the indirect effect and proportion of variance explained by each effect for each mediation model, we will run these two models.
As mentioned above, to test Hypothesis 1, we must first establish that pathogen disgust affects adherence to traditional practices (i.e., CONV_f_z) and that adherence to traditional practices influences religiosity.
# Fitting the models
H1_M1 <- lm(CONV_f_z ~ TDDS_p_z, data = data)
H1_M2 <- lm(CRS_z ~ CONV_f_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H1_M1, which = 1) # For model 1plot(H1_M2, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H1_M1)) # For model 1hist(residuals(H1_M2)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H1_M1))
qqline(residuals(H1_M1), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H1_M2))
qqline(residuals(H1_M2), col = "red")Model 1:
Model 2:
The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$CONV_f_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Traditionalism",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H1_M1, col = "black", lwd = 1)# Summarizing the model
summary(H1_M1)
Call:
lm(formula = CONV_f_z ~ TDDS_p_z, data = data)
Residuals:
Min 1Q Median 3Q Max
-3.4517 -0.6429 0.0402 0.6287 2.8524
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.860e-17 5.893e-02 0.000 1.000
TDDS_p_z 5.260e-03 5.903e-02 0.089 0.929
Residual standard error: 1.002 on 287 degrees of freedom
Multiple R-squared: 2.766e-05, Adjusted R-squared: -0.003457
F-statistic: 0.00794 on 1 and 287 DF, p-value: 0.9291
We can see from the Scatter Plot that there is no relationship between pathogen disgust and traditionalism to speak of (B = .005, t(287) = .089, p = .929), just as with the correlation above.
# Create a scatterplot with the regression line to visualize
plot(data$CONV_f_z, data$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Traditionalism",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H1_M2, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
CRS_z <- as.numeric(data$CRS_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_M2 <- Boot(H1_M2, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_M2) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 5.1338e-17 -6.9336e-04 0.056239 -0.00035888
CONV_f_z 3.1443e-01 5.1009e-05 0.054894 0.31405077
confint(boot_H1_M2) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1074245 0.1127866
CONV_f_z 0.2047469 0.4243772
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, CRS_z)We can see from the Scatter Plot that there is a positive relationship between traditionalism and religiosity.
Based on the estimates, the standardized coefficient for traditionalism is around .31.
The fact that pathogen disgust does not predict traditionalism indicates that traditionalism could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H1_full <- lm(CRS_z ~ TDDS_p_z + CONV_f_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H1_full, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H1_full))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H1_full))
qqline(residuals(H1_full), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H1_full)TDDS_p_z CONV_f_z
1.000028 1.000028
# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_full <- Boot(H1_full, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_full) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 4.1299e-17 -0.00139011 0.054836 -0.0025004
TDDS_p_z 2.4134e-01 -0.00021190 0.057757 0.2419430
CONV_f_z 3.1316e-01 0.00080321 0.050325 0.3138690
confint(boot_H1_full) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1009681 0.1150922
TDDS_p_z 0.1269087 0.3490115
CONV_f_z 0.2120557 0.4121981
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, CRS_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H1_mediate_results <- mediate(model.m = H1_M1,
model.y = H1_full,
treat = "TDDS_p_z",
mediator = "CONV_f_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results for the mediation analysis
summary(H1_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.00165 -0.03517 0.04 0.96
ADE 0.24134 0.13151 0.35 <2e-16 ***
Total Effect 0.24299 0.13360 0.36 <2e-16 ***
Prop. Mediated 0.00678 -0.18265 0.18 0.96
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Based on the bootstrapped CIs not containing zero, both pathogen disgust and traditionalism are significantly positively related to religiosity.
The coefficients for pathogen disgust (.24) and traditionalism (.31) are virtually identical in this model to the models above where they are stand-alone predictors. This suggests that the unique effects of pathogen disgust and traditionalism on religiosity over-and-above each other are completely independent effects.
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through traditionalism (B = .001, BCI = [-.035, .04]).
This result is inconsistent with both the weak and strong version version of Hypothesis 1. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by traditionalism.
As mentioned above, to test Hypothesis 2, we must first establish that pathogen disgust affects out-group avoidance (i.e., GENE_p_z) and that out-group avoidance influences religiosity.
# Fitting the models
H2_M1 <- lm(GENE_p_z ~ TDDS_p_z, data = data)
H2_M2 <- lm(CRS_z ~ GENE_p_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_M1, which = 1) # For model 1plot(H2_M2, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_M1)) # For model 1hist(residuals(H2_M2)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H2_M1))
qqline(residuals(H2_M1), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H2_M2))
qqline(residuals(H2_M2), col = "red")Model 1:
Model 2:
The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 or 1.5 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$GENE_p_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Out-Group Avoidance (GENE_p)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H2_M1, col = "black", lwd = 1)# Summarizing the model
summary(H2_M1)
Call:
lm(formula = GENE_p_z ~ TDDS_p_z, data = data)
Residuals:
Min 1Q Median 3Q Max
-1.9107 -0.8800 0.0804 0.5946 4.0239
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.593e-16 5.892e-02 0.000 1.000
TDDS_p_z 1.440e-02 5.902e-02 0.244 0.807
Residual standard error: 1.002 on 287 degrees of freedom
Multiple R-squared: 0.0002075, Adjusted R-squared: -0.003276
F-statistic: 0.05956 on 1 and 287 DF, p-value: 0.8074
We can see from the Scatter Plot that there is no relationship between pathogen disgust and out-group avoidance to speak of (B = .014, t(287) = .244, p = .807), just as with the correlation above.
# Create a scatterplot with the regression line to visualize
plot(data$GENE_p_z, data$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Out-Group Avoidance (GENE_p)",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H2_M2, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
CRS_z <- as.numeric(data$CRS_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M2 <- Boot(H2_M2, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_M2) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 1.5852e-17 0.00034356 0.059932 0.0020055
GENE_p_z 6.7778e-02 -0.00134646 0.062326 0.0655162
confint(boot_H2_M2) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.11444697 0.1182643
GENE_p_z -0.04812802 0.1936607
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, CRS_z)We can see from the Scatter Plot that there is a positive trend in the relationship between out-group avoidance and religiosity.
Based on the estimates, the standardized coefficient for out-group avoidance is around .068.
The fact that pathogen disgust does not predict out-group avoidance and, in turn, out-group avoidance does not predict religiosity indicates that out-group avoidance could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H2_full <- lm(CRS_z ~ TDDS_p_z + GENE_p_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_full, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_full))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H2_full))
qqline(residuals(H2_full), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H2_full)TDDS_p_z GENE_p_z
1.000208 1.000208
# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_full <- Boot(H2_full, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_full) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 6.4410e-18 -0.00036098 0.058493 -0.0014392
TDDS_p_z 2.4206e-01 -0.00012574 0.059348 0.2420111
GENE_p_z 6.4291e-02 -0.00114233 0.059533 0.0630945
confint(boot_H2_full) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.10737353 0.1269200
TDDS_p_z 0.12066998 0.3534012
GENE_p_z -0.04731359 0.1858159
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, CRS_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H2_mediate_results <- mediate(model.m = H2_M1,
model.y = H2_full,
treat = "TDDS_p_z",
mediator = "GENE_p_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results of the mediation analysis
summary(H2_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.000926 -0.008892 0.01 0.81
ADE 0.242061 0.132686 0.36 <2e-16 ***
Total Effect 0.242987 0.133597 0.36 <2e-16 ***
Prop. Mediated 0.003811 -0.039966 0.06 0.81
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Based on the bootstrapped CIs, pathogen disgust, but not out-group avoidance, is significantly positively related to religiosity.
The coefficients for pathogen disgust (.24) and out-group avoidance (.06) are virtually identical in this model to the models above where they are stand-alone predictors. This suggests that the unique effect of pathogen disgust (and not out-group avoidance) on religiosity over-and-above each other are completely independent.
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through out-group avoidance (B = .001, BCI = [-.009, .01]).
This result is inconsistent with both the weak and strong version version of Hypothesis 2. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by out-group avoidance.
As mentioned above, to test Hypothesis 3, we must first establish that pathogen disgust affects one’s mating strategy (i.e., make it more monogamous; as measured by more restricted sociosexuality; SOI_a_r_z) and that a monogamous mating strategy influences religiosity.
# Fitting the models
H3_M1 <- lm(SOI_a_r_z ~ TDDS_p_z, data = data)
H3_M2 <- lm(CRS_z ~ SOI_a_r_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H3_M1, which = 1) # For model 1plot(H3_M2, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H3_M1)) # For model 1hist(residuals(H3_M2)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H3_M1))
qqline(residuals(H3_M1), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H3_M2))
qqline(residuals(H3_M2), col = "red")Both Models:
The QQ-plots for both models do not indicate that the residuals are normally distributed.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for both models.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$SOI_a_r_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Restricted Sociosexual Attitudes",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H3_M1, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M1 <- Boot(H3_M1, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_M1) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -1.1616e-17 0.0012152 0.059909 0.0021975
TDDS_p_z 2.2087e-01 0.0015469 0.059225 0.2236528
confint(boot_H3_M1) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1256230 0.1076692
TDDS_p_z 0.1039745 0.3299802
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, TDDS_p_z)We can see from the Scatter Plot that there is a positive relationship between pathogen disgust and restricted sociosexual attitudes.
Based on the estimates, the standardized coefficient for pathogen disgust is around .22.
# Create a scatterplot with the regression line to visualize
plot(data$SOI_a_r_z, data$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Restricted Sociosexual Attitudes",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H3_M2, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
CRS_z <- as.numeric(data$CRS_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M2 <- Boot(H3_M2, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_M2) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 2.7721e-17 -0.00099204 0.054590 -0.00413
SOI_a_r_z 4.0934e-01 0.00020098 0.055098 0.41078
confint(boot_H3_M2) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.09552607 0.1229143
SOI_a_r_z 0.29501332 0.5104241
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, CRS_z)We can see from the Scatter Plot that there is a positive trend in the relationship between restricted sociosexual attitudes and religiosity.
Based on the estimates, the standardized coefficient for restricted sociosexual attitudes is around .41.
The fact that pathogen disgust predicts a monogamous mating strategy and, in turn, a monogamous mating strategy predicts religiosity indicates that a monogamous mating strategy could potentially mediate the relationship between pathogen disgust and religiosity. Now we will run the mediation analysis.
# Fitting the full mediation model
H3_full <- lm(CRS_z ~ TDDS_p_z + SOI_a_r_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H3_full, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H3_full))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H3_full))
qqline(residuals(H3_full), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H3_full) TDDS_p_z SOI_a_r_z
1.051286 1.051286
# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_full <- Boot(H3_full, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_full) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 2.1027e-17 -0.00136655 0.054173 -0.0038041
TDDS_p_z 1.6040e-01 -0.00045433 0.056386 0.1585189
SOI_a_r_z 3.7391e-01 0.00013298 0.057102 0.3743473
confint(boot_H3_full) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.09343557 0.1242681
TDDS_p_z 0.05467911 0.2719632
SOI_a_r_z 0.25544697 0.4836007
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, CRS_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H3_mediate_results <- mediate(model.m = H3_M1,
model.y = H3_full,
treat = "TDDS_p_z",
mediator = "SOI_a_r_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results of the mediation analysis
summary(H3_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.0826 0.0384 0.13 <2e-16 ***
ADE 0.1604 0.0563 0.27 0.002 **
Total Effect 0.2430 0.1336 0.36 <2e-16 ***
Prop. Mediated 0.3399 0.1636 0.63 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Mediation Analysis Results:
The ACME indicates that there is a significant indirect effect of pathogen disgust on religiosity through a monogamous mating strategy (B = .08, BCI = [.038, .13]).
The proportion of the total effect mediated (the ratio of the indirect effect to the direct effect) is .34. That is, 34% of the covariance between pathogen disgust and religiosity seems to be able to be accounted for by restricted sociosexual attitudes. Our inferential criterion for concluding full mediation was that this ratio would be at least .80, with the total effect greater than a B = .2. Even the upper bound of the bootstrapped confidence interval for this proportion mediated does not cross this threshold (BCI = [.164, .63]).
We will now re-test Hypothesis 1 while controlling for sex.
# Fitting the models
H1_M1_s <- lm(CONV_f_z ~ TDDS_p_z + sex, data = data)
H1_M2_s <- lm(CRS_z ~ CONV_f_z + sex, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H1_M1_s, which = 1) # For model 1plot(H1_M2_s, which = 1) # For model 2# Extract residuals and fitted values for model 1
residuals_m1 <- residuals(H1_M1_s)
fitted_m1 <- fitted(H1_M1_s)
# Plot residuals vs fitted values for model 1 with different colors for sex
plot(fitted_m1, residuals_m1, col = as.numeric(data$sex),
main = "Residuals vs Fitted Values (Model 1)",
xlab = "Fitted Values", ylab = "Residuals")
abline(h = 0, col = "red")
# Add a legend to the plot
legend("topright", legend = levels(data$sex), col = 1:length(levels(data$sex)), pch = 1)# Drop the residuals and fitted values from the environment
rm(residuals_m1, fitted_m1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H1_M1_s)) # For model 1hist(residuals(H1_M2_s)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H1_M1_s))
qqline(residuals(H1_M1_s), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H1_M2_s))
qqline(residuals(H1_M2_s), col = "red")Model 1:
Model 2:
We will conduct a bootstrapped tests of the predictors for both models.
Multicollinearity: Variance inflation factor (VIF) for both models
# Calculating the VIF for both models
vif(H1_M1_s)TDDS_p_z sex
1.038981 1.038981
vif(H1_M2_s)CONV_f_z sex
1.004388 1.004388
# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_M1_s <- Boot(H1_M1_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_M1_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 0.068890 9.5639e-04 0.083729 0.071961
TDDS_p_z 0.018766 -8.7770e-04 0.067435 0.017903
sexFemale -0.139225 3.8768e-06 0.118828 -0.140871
confint(boot_H1_M1_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1076349 0.22374169
TDDS_p_z -0.1137392 0.15325593
sexFemale -0.3708032 0.09485534
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, TDDS_p_z, sex)Based on the estimate and BCI for pathogen disgust (B = .018, BCI = [-.114, .153]), controlling for sex does not change whether there is a positive relationship between pathogen disgust and traditionalism.
# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
CRS_z <- as.numeric(data$CRS_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_M2_s <- Boot(H1_M2_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_M2_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -0.18812 -0.00399574 0.074103 -0.19341
CONV_f_z 0.32702 0.00040075 0.053030 0.32711
sexFemale 0.38019 0.00631972 0.113040 0.38812
confint(boot_H1_M2_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.3170446 -0.02384113
CONV_f_z 0.2199292 0.42953209
sexFemale 0.1614658 0.60697545
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, CRS_z, sex)The fact that pathogen disgust does not predict traditionalism indicates that traditionalism could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H1_full_s <- lm(CRS_z ~ TDDS_p_z + CONV_f_z + sex, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H1_full_s, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H1_full_s))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H1_full_s))
qqline(residuals(H1_full_s), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a tendency for higher frequencies of residuals at just around 1 SD.
We will conduct bootstrapping for confidence intervals to make inferences.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H1_full_s)TDDS_p_z CONV_f_z sex
1.039334 1.004730 1.043866
# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_full_s <- Boot(H1_full_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_full_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -0.14722 -0.00473112 0.074568 -0.15385
TDDS_p_z 0.21242 -0.00125691 0.057959 0.21105
CONV_f_z 0.32316 0.00096636 0.049858 0.32491
sexFemale 0.29753 0.00669264 0.110568 0.30430
confint(boot_H1_full_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.27878504 0.02449477
TDDS_p_z 0.09861968 0.32435475
CONV_f_z 0.21945136 0.41880569
sexFemale 0.07606250 0.50767084
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, CRS_z, TDDS_p_z, sex)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H1_s_mediate_results <- mediate(model.m = H1_M1_s,
model.y = H1_full_s,
treat = "TDDS_p_z",
mediator = "CONV_f_z",
covariates = "sex",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results for the mediation analysis
summary(H1_s_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
(Inference Conditional on the Covariate Values Specified in `covariates')
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.00606 -0.03225 0.05 0.76
ADE 0.21242 0.10577 0.33 <2e-16 ***
Total Effect 0.21849 0.10907 0.35 <2e-16 ***
Prop. Mediated 0.02776 -0.19673 0.23 0.76
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through traditionalism (B = .006, BCI = [-.032, .05]).
This result is inconsistent with both the weak and strong version version of Hypothesis 1. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by traditionalism, even when controlling for sex.
We will now re-test Hypothesis 2 while controlling for sex.
# Fitting the models
H2_M1_s <- lm(GENE_p_z ~ TDDS_p_z + sex, data = data)
H2_M2_s <- lm(CRS_z ~ GENE_p_z + sex, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_M1_s, which = 1) # For model 1plot(H2_M2_s, which = 1) # For model 2# Extract residuals and fitted values for model 1
residuals_m1 <- residuals(H2_M1_s)
fitted_m1 <- fitted(H2_M1_s)
# Plot residuals vs fitted values for model 1 with different colors for sex
plot(fitted_m1, residuals_m1, col = as.numeric(data$sex),
main = "Residuals vs Fitted Values (Model 1)",
xlab = "Fitted Values", ylab = "Residuals")
abline(h = 0, col = "red")
# Add a legend to the plot
legend("topright", legend = levels(data$sex), col = 1:length(levels(data$sex)), pch = 1)# Drop the residuals and fitted values from the environment
rm(residuals_m1, fitted_m1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_M1_s)) # For model 1hist(residuals(H2_M2_s)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H2_M1_s))
qqline(residuals(H2_M1_s), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H2_M2_s))
qqline(residuals(H2_M2_s), col = "red")Model 1:
Model 2:
We will conduct a bootstrapped tests of the predictors for both models.
Multicollinearity: Variance inflation factor (VIF) for both models
# Calculating the VIF for both models
vif(H2_M1_s)TDDS_p_z sex
1.038981 1.038981
vif(H2_M2_s)GENE_p_z sex
1.013159 1.013159
# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M1_s <- Boot(H2_M1_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_M1_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 0.119846 -0.00150248 0.087176 0.118491
TDDS_p_z 0.037901 -0.00194585 0.061237 0.036581
sexFemale -0.242207 -0.00069038 0.120490 -0.243567
confint(boot_H2_M1_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.04583828 0.29032402
TDDS_p_z -0.08378796 0.15652933
sexFemale -0.47809343 0.01058525
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, TDDS_p_z, sex)Based on the estimate and BCI for pathogen disgust (B = .038, BCI = [-.084, .157]), while controlling for sex there is not a relationship between pathogen disgust and out-group avoidance.
# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
CRS_z <- as.numeric(data$CRS_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M2_s <- Boot(H2_M2_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_M2_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -0.17669 -0.00266419 0.078782 -0.183267
GENE_p_z 0.08816 -0.00091022 0.061117 0.085621
sexFemale 0.35710 0.00584804 0.118083 0.362284
confint(boot_H2_M2_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.30605613 0.0008240374
GENE_p_z -0.02455483 0.2142069266
sexFemale 0.13007365 0.5932166829
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, CRS_z, sex)The fact that pathogen disgust does not predict out-group avoidance and that out-group avoidance does not predict religiosity indicates that traditionalism could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H2_full_s <- lm(CRS_z ~ TDDS_p_z + GENE_p_z + sex, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_full_s, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_full_s))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H2_full_s))
qqline(residuals(H2_full_s), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a tendency for higher frequencies of residuals at just above -1 SD and just above 1 SD.
We will conduct bootstrapping for confidence intervals to make inferences.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H2_full_s)TDDS_p_z GENE_p_z sex
1.040438 1.014580 1.053911
# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_full_s <- Boot(H2_full_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_full_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -0.134569 -0.0036411 0.079129 -0.142816
TDDS_p_z 0.215448 -0.0012646 0.060393 0.214011
GENE_p_z 0.080198 -0.0005672 0.058824 0.078091
sexFemale 0.271962 0.0066281 0.116184 0.278642
confint(boot_H2_full_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.26815673 0.03305607
TDDS_p_z 0.09131601 0.33257881
GENE_p_z -0.02763633 0.20025791
sexFemale 0.03676185 0.49832392
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, CRS_z, TDDS_p_z, sex)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H2_s_mediate_results <- mediate(model.m = H2_M1_s,
model.y = H2_full_s,
treat = "TDDS_p_z",
mediator = "GENE_p_z",
covariates = "sex",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results for the mediation analysis
summary(H2_s_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
(Inference Conditional on the Covariate Values Specified in `covariates')
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.00304 -0.00715 0.02 0.58
ADE 0.21545 0.10555 0.34 <2e-16 ***
Total Effect 0.21849 0.10907 0.35 <2e-16 ***
Prop. Mediated 0.01391 -0.03591 0.09 0.58
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through out-group avoidance (B = .003, BCI = [-.007, .02]).
This result is inconsistent with both the weak and strong version version of Hypothesis 2. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by out-group avoidance, even when controlling for sex.
We will now re-test Hypothesis 3 while controlling for sex.
# Fitting the models
H3_M1_s <- lm(SOI_a_r_z ~ TDDS_p_z + sex, data = data)
H3_M2_s <- lm(CRS_z ~ SOI_a_r_z + sex, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H3_M1_s, which = 1) # For model 1plot(H3_M2_s, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H3_M1_s)) # For model 1hist(residuals(H3_M2_s)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H3_M1_s))
qqline(residuals(H3_M1_s), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H3_M2_s))
qqline(residuals(H3_M2_s), col = "red")Both Models:
We will conduct a bootstrapped tests of the predictors for both models.
Multicollinearity: Variance inflation factor (VIF) for both models
# Calculating the VIF for both models
vif(H3_M1_s)TDDS_p_z sex
1.038981 1.038981
vif(H3_M2_s)SOI_a_r_z sex
1.038623 1.038623
# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M1_s <- Boot(H3_M1_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_M1_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -0.15403 -0.00051909 0.083356 -0.15251
TDDS_p_z 0.19067 0.00065110 0.059984 0.19025
sexFemale 0.31129 0.00377787 0.113956 0.31817
confint(boot_H3_M1_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.3198589 0.007813617
TDDS_p_z 0.0761788 0.308853893
sexFemale 0.0794630 0.540941360
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, TDDS_p_z, sex)# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
CRS_z <- as.numeric(data$CRS_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M2_s <- Boot(H3_M2_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_M2_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -0.092212 -0.00303183 0.078827 -0.10045
SOI_a_r_z 0.391337 -0.00026862 0.055955 0.39226
sexFemale 0.186359 0.00421532 0.113045 0.18804
confint(boot_H3_M2_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.23088030 0.08466747
SOI_a_r_z 0.27130629 0.49262933
sexFemale -0.02599912 0.41658472
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, CRS_z, sex)# Fitting the full mediation model
H3_full_s <- lm(CRS_z ~ TDDS_p_z + SOI_a_r_z + sex, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H3_full_s, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H3_full_s))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H3_full_s))
qqline(residuals(H3_full_s), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not quite normally distributed.
We will conduct bootstrapping for confidence intervals to make inferences.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H3_full_s) TDDS_p_z SOI_a_r_z sex
1.078165 1.077794 1.065178
# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
sex <- data$sex
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_full_s <- Boot(H3_full_s, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_full_s) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -0.069064 -0.0038943 0.078950 -0.078932
TDDS_p_z 0.149296 -0.0010695 0.057400 0.148526
SOI_a_r_z 0.362880 -0.0001937 0.057988 0.361283
sexFemale 0.139577 0.0051513 0.113132 0.146277
confint(boot_H3_full_s) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.20809145 0.1084253
TDDS_p_z 0.03664188 0.2590547
SOI_a_r_z 0.24555603 0.4757293
sexFemale -0.08302237 0.3537507
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, CRS_z, TDDS_p_z, sex)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H3_s_mediate_results <- mediate(model.m = H3_M1_s,
model.y = H3_full_s,
treat = "TDDS_p_z",
mediator = "SOI_a_r_z",
covariates = "sex",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results for the mediation analysis
summary(H3_s_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
(Inference Conditional on the Covariate Values Specified in `covariates')
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.0692 0.0266 0.12 <2e-16 ***
ADE 0.1493 0.0431 0.27 0.008 **
Total Effect 0.2185 0.1091 0.35 <2e-16 ***
Prop. Mediated 0.3167 0.1362 0.62 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Mediation Analysis Results:
The ACME indicates that there is a significant indirect effect of pathogen disgust on religiosity through a monogamous mating strategy (B = .069, BCI = [.027, .12]).
The proportion of the total effect mediated (the ratio of the indirect effect to the direct effect) is .32. That is, 32% of the covariance between pathogen disgust and religiosity seems to be able to be accounted for by restricted sociosexual attitudes. Our inferential criterion for concluding full mediation was that this ratio would be at least .80, with the total effect greater than a B = .2. Even the upper bound of the bootstrapped confidence interval for this proportion mediated does not cross this threshold (BCI = [.136, .62]).
Now we will retest each hypothesis with religious importance as the outcome variable, as a check for robustness.
As mentioned above, to test Hypothesis 1, we must first establish that pathogen disgust affects adherence to traditional practices (i.e., CONV_f_z) and that adherence to traditional practices influences religious importance.
# Fitting the models
H1_M1_RI <- lm(CONV_f_z ~ TDDS_p_z, data = data)
H1_M2_RI <- lm(religious_importance_z ~ CONV_f_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H1_M1_RI, which = 1) # For model 1plot(H1_M2_RI, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H1_M1_RI)) # For model 1hist(residuals(H1_M2_RI)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H1_M1_RI))
qqline(residuals(H1_M1_RI), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H1_M2_RI))
qqline(residuals(H1_M2_RI), col = "red")Model 1:
Model 2:
The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$CONV_f_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Traditionalism",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H1_M1_RI, col = "black", lwd = 1)# Summarizing the model
summary(H1_M1_RI)
Call:
lm(formula = CONV_f_z ~ TDDS_p_z, data = data)
Residuals:
Min 1Q Median 3Q Max
-3.4517 -0.6429 0.0402 0.6287 2.8524
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7.860e-17 5.893e-02 0.000 1.000
TDDS_p_z 5.260e-03 5.903e-02 0.089 0.929
Residual standard error: 1.002 on 287 degrees of freedom
Multiple R-squared: 2.766e-05, Adjusted R-squared: -0.003457
F-statistic: 0.00794 on 1 and 287 DF, p-value: 0.9291
We can see from the Scatter Plot that there is no relationship between pathogen disgust and traditionalism to speak of (B = .005, t(287) = .089, p = .929), just as with the correlation above.
# Create a scatterplot with the regression line to visualize
plot(data$CONV_f_z, data$religious_importance_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Traditionalism",
ylab = "Religious Importance",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H1_M2_RI, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
religious_importance_z <- as.numeric(data$religious_importance_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_M2_RI <- Boot(H1_M2_RI, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_M2_RI) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 2.3298e-17 4.6614e-06 0.055910 -0.0027938
CONV_f_z 3.3199e-01 1.4028e-03 0.055438 0.3338146
confint(boot_H1_M2_RI) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1002662 0.1200648
CONV_f_z 0.2204035 0.4370547
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, religious_importance_z)We can see from the Scatter Plot that there is a positive relationship between traditionalism and religious importance.
Based on the estimates, the standardized coefficient for traditionalism is around .33.
The fact that pathogen disgust does not predict traditionalism indicates that traditionalism could not mediate the relationship between pathogen disgust and religious importance, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H1_full_RI <- lm(religious_importance_z ~ TDDS_p_z + CONV_f_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H1_full_RI, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H1_full_RI))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H1_full_RI))
qqline(residuals(H1_full_RI), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just below 1 SD.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H1_full_RI)TDDS_p_z CONV_f_z
1.000028 1.000028
# Adding the variables to the broader environment so the Boot function will run
CONV_f_z <- as.numeric(data$CONV_f_z)
religious_importance_z <- as.numeric(data$religious_importance_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H1_full_RI <- Boot(H1_full_RI, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H1_full_RI) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 1.4280e-17 -0.00064525 0.054844 -0.0038586
TDDS_p_z 2.1682e-01 0.00030940 0.057438 0.2166807
CONV_f_z 3.3085e-01 0.00205738 0.051646 0.3331901
confint(boot_H1_full_RI) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.09907331 0.1172830
TDDS_p_z 0.10874782 0.3291841
CONV_f_z 0.22029722 0.4256569
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(CONV_f_z, religious_importance_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H1_RI_mediate_results <- mediate(model.m = H1_M1_RI,
model.y = H1_full_RI,
treat = "TDDS_p_z",
mediator = "CONV_f_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results for the mediation analysis
summary(H1_RI_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.00174 -0.03773 0.04 0.96
ADE 0.21682 0.10884 0.33 <2e-16 ***
Total Effect 0.21856 0.11198 0.34 <2e-16 ***
Prop. Mediated 0.00796 -0.22993 0.21 0.96
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Based on the bootstrapped CIs not containing zero, both pathogen disgust and traditionalism are significantly positively related to religious importance.
The coefficients for pathogen disgust (.22) and traditionalism (.33) are similar to the model with the CRS.
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religious importance through traditionalism (B = .002, BCI = [-.038, .04]).
This result is inconsistent with both the weak and strong version version of Hypothesis 1. That is, there is no evidence here that the relationship between pathogen disgust and religious importance is mediated by traditionalism.
As mentioned above, to test Hypothesis 2, we must first establish that pathogen disgust affects out-group avoidance (i.e., GENE_p_z) and that out-group avoidance influences religious importance.
# Fitting the models
H2_M1_RI <- lm(GENE_p_z ~ TDDS_p_z, data = data)
H2_M2_RI <- lm(religious_importance_z ~ GENE_p_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_M1_RI, which = 1) # For model 1plot(H2_M2_RI, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_M1_RI)) # For model 1hist(residuals(H2_M2_RI)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H2_M1_RI))
qqline(residuals(H2_M1_RI), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H2_M2_RI))
qqline(residuals(H2_M2_RI), col = "red")Model 1:
Model 2:
The histogram and QQ-plot of the residuals of Model 2 indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 or 1.5 SD.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for Model 2.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$GENE_p_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Out-Group Avoidance (GENE_p)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H2_M1_RI, col = "black", lwd = 1)# Summarizing the model
summary(H2_M1_RI)
Call:
lm(formula = GENE_p_z ~ TDDS_p_z, data = data)
Residuals:
Min 1Q Median 3Q Max
-1.9107 -0.8800 0.0804 0.5946 4.0239
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.593e-16 5.892e-02 0.000 1.000
TDDS_p_z 1.440e-02 5.902e-02 0.244 0.807
Residual standard error: 1.002 on 287 degrees of freedom
Multiple R-squared: 0.0002075, Adjusted R-squared: -0.003276
F-statistic: 0.05956 on 1 and 287 DF, p-value: 0.8074
We can see from the Scatter Plot that there is no relationship between pathogen disgust and out-group avoidance to speak of (B = .014, t(287) = .244, p = .807)—same as the model above.
# Create a scatterplot with the regression line to visualize
plot(data$GENE_p_z, data$religious_importance_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Out-Group Avoidance (GENE_p)",
ylab = "Religious Importance",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H2_M2_RI, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
religious_importance_z <- as.numeric(data$religious_importance_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M2_RI <- Boot(H2_M2_RI, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_M2_RI) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -1.4709e-17 0.00069972 0.059421 0.00035152
GENE_p_z 7.4934e-02 -0.00254865 0.062068 0.07214367
confint(boot_H2_M2_RI) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.11600555 0.1186471
GENE_p_z -0.04012943 0.1972286
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, religious_importance_z)We can see from the Scatter Plot that there is a positive trend in the relationship between out-group avoidance and religious importance.
Based on the estimates, the standardized coefficient for traditionalism is around .07.
The fact that pathogen disgust does not predict out-group avoidance and, in turn, out-group avoidance does not predict religious importance indicates that out-group avoidance could not mediate the relationship between pathogen disgust and religious importance, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H2_full_RI <- lm(religious_importance_z ~ TDDS_p_z + GENE_p_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_full_RI, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_full_RI))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H2_full_RI))
qqline(residuals(H2_full_RI), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H2_full_RI)TDDS_p_z GENE_p_z
1.000208 1.000208
# Adding the variables to the broader environment so the Boot function will run
GENE_p_z <- as.numeric(data$GENE_p_z)
religious_importance_z <- as.numeric(data$religious_importance_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_full_RI <- Boot(H2_full_RI, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_full_RI) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -2.3166e-17 5.6297e-05 0.058364 5.6535e-05
TDDS_p_z 2.1752e-01 2.5512e-04 0.059996 2.1741e-01
GENE_p_z 7.1801e-02 -2.2410e-03 0.059293 6.9262e-02
confint(boot_H2_full_RI) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.11278531 0.1217364
TDDS_p_z 0.09532824 0.3302613
GENE_p_z -0.03957865 0.1880959
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_z, religious_importance_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H2_RI_mediate_results <- mediate(model.m = H2_M1_RI,
model.y = H2_full_RI,
treat = "TDDS_p_z",
mediator = "GENE_p_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results of the mediation analysis
summary(H2_RI_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.00103 -0.00902 0.01 0.83
ADE 0.21752 0.11006 0.34 <2e-16 ***
Total Effect 0.21856 0.11198 0.34 <2e-16 ***
Prop. Mediated 0.00473 -0.04487 0.07 0.83
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Based on the bootstrapped CIs, pathogen disgust, but not out-group avoidance, is significantly positively related to religious importance.
The coefficients for pathogen disgust (.22) and out-group avoidance (.07) are virtually identical in the model with CRS.
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religious importance through out-group avoidance (B = .001, BCI = [-.009, .01]).
This result is inconsistent with both the weak and strong version version of Hypothesis 2. That is, there is no evidence here that the relationship between pathogen disgust and religious importance is mediated by out-group avoidance.
As mentioned above, to test Hypothesis 3, we must first establish that pathogen disgust affects one’s mating strategy (i.e., make it more monogamous; as measured by more restricted sociosexuality; SOI_a_r_z) and that a monogamous mating strategy influences religious importance.
# Fitting the models
H3_M1_RI <- lm(SOI_a_r_z ~ TDDS_p_z, data = data)
H3_M2_RI <- lm(religious_importance_z ~ SOI_a_r_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H3_M1_RI, which = 1) # For model 1plot(H3_M2_RI, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H3_M1_RI)) # For model 1hist(residuals(H3_M2_RI)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H3_M1_RI))
qqline(residuals(H3_M1_RI), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H3_M2_RI))
qqline(residuals(H3_M2_RI), col = "red")Both Models:
The QQ-plots for both models do not indicate that the residuals are normally distributed.
Because of this failure of this assumption, we will conduct a bootstrapped test of the predictor for both models.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$SOI_a_r_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Restricted Sociosexual Attitudes",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H3_M1_RI, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M1_RI <- Boot(H3_M1_RI, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_M1_RI) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -1.1616e-17 0.0012152 0.059909 0.0021975
TDDS_p_z 2.2087e-01 0.0015469 0.059225 0.2236528
confint(boot_H3_M1_RI) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1256230 0.1076692
TDDS_p_z 0.1039745 0.3299802
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, TDDS_p_z)We can see from the Scatter Plot that there is a positive relationship between pathogen disgust and restricted sociosexual attitudes.
Based on the estimates, the standardized coefficient for pathogen disgust is around .22.
# Create a scatterplot with the regression line to visualize
plot(data$SOI_a_r_z, data$religious_importance_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Restricted Sociosexual Attitudes",
ylab = "Religious Importance",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H3_M2_RI, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
religious_importance_z <- as.numeric(data$religious_importance_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_M2_RI <- Boot(H3_M2_RI, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_M2_RI) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -1.7538e-18 -0.00046013 0.054763 -0.00017089
SOI_a_r_z 3.8592e-01 -0.00012983 0.056140 0.38640935
confint(boot_H3_M2_RI) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1039343 0.1170591
SOI_a_r_z 0.2740123 0.4924135
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, religious_importance_z)We can see from the Scatter Plot that there is a positive trend in the relationship between restricted sociosexual attitudes and religious importance.
Based on the estimates, the standardized coefficient for restricted sociosexual attitudes is around .39.
The fact that pathogen disgust predicts a monogamous mating strategy and, in turn, a monogamous mating strategy predicts religious importance indicates that a monogamous mating strategy could potentially mediate the relationship between pathogen disgust and religious importance. Now we will run the mediation analysis.
# Fitting the full mediation model
H3_full_RI <- lm(religious_importance_z ~ TDDS_p_z + SOI_a_r_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H3_full_RI, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H3_full_RI))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H3_full_RI))
qqline(residuals(H3_full_RI), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H3_full_RI) TDDS_p_z SOI_a_r_z
1.051286 1.051286
# Adding the variables to the broader environment so the Boot function will run
SOI_a_r_z <- as.numeric(data$SOI_a_r_z)
religious_importance_z <- as.numeric(data$religious_importance_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H3_full_RI <- Boot(H3_full_RI, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H3_full_RI) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) -7.6035e-18 -0.00082361 0.054498 -0.0035079
TDDS_p_z 1.4015e-01 0.00033273 0.055713 0.1394299
SOI_a_r_z 3.5496e-01 -0.00063970 0.057809 0.3543830
confint(boot_H3_full_RI) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.09587216 0.1273944
TDDS_p_z 0.03656495 0.2534092
SOI_a_r_z 0.23746177 0.4667500
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(SOI_a_r_z, religious_importance_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H3_RI_mediate_results <- mediate(model.m = H3_M1_RI,
model.y = H3_full_RI,
treat = "TDDS_p_z",
mediator = "SOI_a_r_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results of the mediation analysis
summary(H3_RI_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME 0.0784 0.0368 0.13 <2e-16 ***
ADE 0.1402 0.0376 0.25 0.002 **
Total Effect 0.2186 0.1120 0.34 <2e-16 ***
Prop. Mediated 0.3587 0.1757 0.71 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Mediation Analysis Results:
The ACME indicates that there is a significant indirect effect of pathogen disgust on religious importance through a monogamous mating strategy (B = .078, BCI = [.037, .13]).
The proportion of the total effect mediated (the ratio of the indirect effect to the direct effect) is .36. That is, 36% of the covariance between pathogen disgust and religious importance seems to be able to be accounted for by restricted sociosexual attitudes. Our inferential criterion for concluding full mediation was that this ratio would be at least .80, with the total effect greater than a B = .2. Even the upper bound of the bootstrapped confidence interval for this proportion mediated does not cross this threshold (BCI = [.176, .71]).
In the Internal Reliability analysis of the Preferences Subscale of SFGENE-7, the drop alpha statistics indicated that removing Item 1 from the subscale would improve the internal reliability from an alpha of around .5 to around .8. Given that a lack of internal reliability may attenuate relationships between variables in correlational analyses, I will now conduct the same analysis for Hypothesis 2 but with Item 1 removed.
As mentioned above, to test Hypothesis 2, we must first establish that pathogen disgust affects out-group avoidance (i.e., GENE_p_2_z) and that out-group avoidance influences religiosity.
# Fitting the models
H2_M1_2 <- lm(GENE_p_2_z ~ TDDS_p_z, data = data)
H2_M2_2 <- lm(CRS_z ~ GENE_p_2_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_M1_2, which = 1) # For model 1plot(H2_M2_2, which = 1) # For model 2Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_M1_2)) # For model 1hist(residuals(H2_M2_2)) # For model 2# QQ-plot for normality of residuals of Model 1
qqnorm(residuals(H2_M1_2))
qqline(residuals(H2_M1_2), col = "red")# QQ-plot for normality of residuals of Model 2
qqnorm(residuals(H2_M2_2))
qqline(residuals(H2_M2_2), col = "red")Both Models:
The histograms and QQ-plots indicate that for both models the residuals are not normally distributed, although in different ways.
We will need to bootstrapping for testing the predictors for these models.
Multicollinearity: Not applicable because each model only has one predictor.
# Create a scatterplot with the regression line to visualize
plot(data$TDDS_p_z, data$GENE_p_2_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Pathogen Disgust",
ylab = "Out-Group Avoidance (GENE_p without Item 1)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H2_M1_2, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
GENE_p_2_z <- as.numeric(data$GENE_p_2_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M1_2 <- Boot(H2_M1_2, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_M1_2) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 7.0301e-16 -0.0014508 0.058715 -0.0040643
TDDS_p_z -2.5360e-02 -0.0020741 0.059085 -0.0284347
confint(boot_H2_M1_2) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1071551 0.1267461
TDDS_p_z -0.1333744 0.1058820
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_2_z, TDDS_p_z)We can see from the Scatter Plot that there is no relationship between pathogen disgust and out-group avoidance to speak of (B = -.025, BCI = [-.133, .106])—just as with the analysis above.
# Create a scatterplot with the regression line to visualize
plot(data$GENE_p_2_z, data$CRS_z,
main = "Scatter Plot with Linear Regression Line",
xlab = "Out-Group Avoidance (GENE_p without Item 1)",
ylab = "Religiosity (CRS)",
pch = 16, # For solid circle points
col = "black")
# Add regression line
abline(H2_M2_2, col = "black", lwd = 1)# Adding the variables to the broader environment so the Boot function will run
GENE_p_2_z <- as.numeric(data$GENE_p_2_z)
CRS_z <- as.numeric(data$CRS_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_M2_2 <- Boot(H2_M2_2, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_M2_2) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 2.6680e-17 -0.00011423 0.060017 0.0015515
GENE_p_2_z 1.3388e-05 -0.00097717 0.064575 -0.0015445
confint(boot_H2_M2_2) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1135425 0.1163811
GENE_p_2_z -0.1190670 0.1228603
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_2_z, CRS_z)We can see from the Scatter Plot that there is a positive trend in the relationship between out-group avoidance and religiosity.
Based on the estimates, the standardized coefficient for out-group avoidance is around B = .00001 (BCI = [-.119, .123]).
The fact that pathogen disgust does not predict out-group avoidance and, in turn, out-group avoidance does not predict religiosity indicates that out-group avoidance could not mediate the relationship between pathogen disgust and religiosity, but I will run the analysis that we planned anyway just to have the results.
# Fitting the full mediation model
H2_2_full <- lm(CRS_z ~ TDDS_p_z + GENE_p_2_z, data = data)Linearity and Heteroscedasticity: Assessed visually through plotting residual vs. fitted values
# Plotting the residual vs. fitted values
plot(H2_2_full, which = 1)Normality: Histogram and QQ-plot (because Shapiro-Wilk test would be overpowered)
# Histogram for normality of residuals
hist(residuals(H2_2_full))# QQ-plot for normality of residuals of the full model
qqnorm(residuals(H2_2_full))
qqline(residuals(H2_2_full), col = "red")The histogram and QQ-plot of the residuals of the model indicate that the residuals are not normally distributed. That is, there is a strong tendency for higher frequencies of residuals at both just above -1 SD and just above 1 SD. The two clusters of similar errors in prediction may be due to sex differences in religiosity. In exploratory analyses, we will add sex as a control variable.
Multicollinearity: Variance Inflation Factor (VIF) from car package
# Calculating VIF
vif(H2_2_full) TDDS_p_z GENE_p_2_z
1.000644 1.000644
# Adding the variables to the broader environment so the Boot function will run
GENE_p_2_z <- as.numeric(data$GENE_p_2_z)
CRS_z <- as.numeric(data$CRS_z)
TDDS_p_z <- as.numeric(data$TDDS_p_z)
# Setting the seed
set.seed(1042557)
# Getting the bootstrapped estimates and CIs (using car package rather than boot for simplicity)
boot_H2_2_full <- Boot(H2_2_full, R = 1000) # 1,000 resamples
# Summarizing the estimates and CIs
summary(boot_H2_2_full) # Getting the coefficients
Number of bootstrap replications R = 1000
original bootBias bootSE bootMed
(Intercept) 1.2339e-17 -0.00077435 0.058551 -0.0013948
TDDS_p_z 2.4314e-01 -0.00072365 0.059883 0.2433389
GENE_p_2_z 6.1794e-03 -0.00060177 0.060073 0.0069054
confint(boot_H2_2_full) # Getting the BCIsBootstrap bca confidence intervals
2.5 % 97.5 %
(Intercept) -0.1079394 0.1232773
TDDS_p_z 0.1198619 0.3565631
GENE_p_2_z -0.1149113 0.1168755
# Removing the vectors in the environment that the Boot package required outside of the data.frame
rm(GENE_p_2_z, CRS_z, TDDS_p_z)
# Using the mediate function from the mediation package to estimate the indirect effect with 95% bootstrapped confidence intervals
H2_2_mediate_results <- mediate(model.m = H2_M1_2,
model.y = H2_2_full,
treat = "TDDS_p_z",
mediator = "GENE_p_2_z",
boot = TRUE, # Bootstrap the CI for the indirect effect
sims = 1000) # Do 1000 resamplesRunning nonparametric bootstrap
# Printing the results of the mediation analysis
summary(H2_2_mediate_results)
Causal Mediation Analysis
Nonparametric Bootstrap Confidence Intervals with the Percentile Method
Estimate 95% CI Lower 95% CI Upper p-value
ACME -0.000157 -0.006528 0.01 0.99
ADE 0.243143 0.131919 0.36 <2e-16 ***
Total Effect 0.242987 0.133597 0.36 <2e-16 ***
Prop. Mediated -0.000645 -0.031489 0.04 0.99
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sample Size Used: 289
Simulations: 1000
Multiple Regression Results:
Based on the bootstrapped CIs, pathogen disgust, but not out-group avoidance, is significantly positively related to religiosity.
The coefficients for pathogen disgust (.24) and out-group avoidance (.006) are similar to their values with Item 1 included.
Mediation Analysis Results:
The ACME indicates that there is no significant indirect effect of pathogen disgust on religiosity through out-group avoidance (B = -.0002, BCI = [-.006, .01]).
This result is inconsistent with both the weak and strong version version of Hypothesis 2. That is, there is no evidence here that the relationship between pathogen disgust and religiosity is mediated by out-group avoidance.
Now we will save the final version of the data frame (data) as a data/final_data.csv.
# Saving the data
write.csv(data, "./data/final_data.csv")